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The predictive power of the ab initio Bethe-Salpeter equation (BSE) approach, rigorously based on many- 
body Green’s function theory but incorporating information from density functional theory, has already 
been demonstrated for the optical gaps and spectra of solid-state systems. Interest in photoactive hybrid 
organic/inorganic systems has recently increased, and so has the use of the BSE for computing neutral 
excitations of organic molecules. However, no systematic benchmarks of the BSE for neutral electronic 
excitations of organic molecules exist. Here, we study the performance of the BSE for the 28 small molecules 
in Thiel’s widely-used time-dependent density functional theory benchmark set [M. Schreiber et al J. Chem. 
Phys. 128 , 134110 (2008)]. We observe that the BSE produces results that depend critically on the mean- 
field starting point employed in the perturbative approach. We find that this starting point dependence 
is mainly introduced through the quasiparticle energies obtained at the intermediate GW step, and that 
with a judicious choice of starting mean-field, singlet excitation energies obtained from BSE are in excellent 
quantitative agreement with higher-level wavefunction methods. The quality of the triplet excitations is 
slightly less satisfactory. 


I. INTRODUCTION 

Optical properties are of broad fundamental and prac¬ 
tical interest. Eor example, they determine the color of 
everyday objects; they dictate the absorption and trans¬ 
fer of photons by and between chromophores embedded 
in protein environments; and they control the fluores¬ 
cence behavior of molecules used as markers in biomedi¬ 
cal imaging applications. These diverse phenomena (and 
many others) are in fact united by the same underlying 
quantum mechanics describing electronic excitations and 
their consequences. 

The ability to reliably and quantitatively predict these 
neutral excitations from computer calculations is an im¬ 
portant goal, and several competing methods for this pur¬ 
pose are in use today. Besides the wavefunction-based 
methods which are computationally-intensive and thus 
often limited to relatively small systems, two main for¬ 
malisms are present in the literature: time-dependent 
density functional theory (TD-DET)^’^ and the Bethe- 
Salpeter equation (BSE)^ approach from many-body per¬ 
turbation theory. Both approaches have been widely used 
but for different classes of systems: TD-DET is primarily 
for molecules, and BSE for solids.^ 

TD-DET has thus far been applied with great success 
to the calculations of the low-lying excitations of isolated 
molecules. It is a computationally attractive method that 
can be used efficiently even for relatively large systems, 
and performs particularly well when paired with a hy¬ 
brid exchange correlation (xc) functional.^ However, TD- 


DET, with standard hybrids, can be challenged by Ryd¬ 
berg final states and charge transfer excitations.^’^ These 
failures can be tempered by employing xc approxima¬ 
tions that contain the full long-range exact exchange con¬ 
tribution, as do the tuned-CAM-BL3YP,^ BNL,^’^^ and 
OT-RSH functionals.^^ However, and despite interesting 
proposals^^”^^, for solids no purely ab initio approxima¬ 
tion to TD-DET is available to date, since the content of 
long-range exact exchange in solids should be modulated 
by the system-dependent dielectric constant. 

Likewise the BSE, whose solution is the two-particle 
electron-hole Green’s function, also accounts for the 
screened electron-hole interaction; the BSE framework 
has been shown to be extremely successful in predict¬ 
ing the optical spectra of bulk solids^^"^^ and of low¬ 
dimensional materials.^^ In part inspired by recent in¬ 
terest in organic-based energy conversion materials, BSE 
has begun to be applied to finite organic molecular sys¬ 
tems as well.^^^^^ Relative to standard contemporary 
TD-DET approaches, the BSE method has many attrac¬ 
tive features: through the ab initio calculation of the 
screened Goulomb interaction, the electron-hole interac¬ 
tion has the correct asymptotic behavior independent of 
the system, be it a bulk solid, a low-dimensional nanos¬ 
tructure or polymer, or a molecule. This feature results, 
for instance, in a correct description of charge transfer 
excitations in molecules. 

Additionally, the description of neutral excitations 
within the BSE are built upon a foundation of 
charged excitation energies, corresponding to elec- 
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tron addition or removal, determined within the GW 
approximation^’^^’^^. The GW approach is known to 
yield much more accurate values of the fundamental 
(or quasiparticle) gap energy for a variety of systems 
than, e.g., DFT. This is in contrast with TD-DFT 
for which underlying Kohn-Sham eigenvalues have little 
physical meaning.Only the highest occupied molecu¬ 
lar orbital (HOMO) and the lowest unoccupied molecu¬ 
lar orbital (LUMO) energies can be safely interpreted as 
the negative of ionization potential (IP) and the neg¬ 
ative of electron affinity (EA) in a generalized Kohn- 
Sham scheme (gKS) (including hybrid functionals). 

All the other eigenvalues are not, strictly speaking, 
observables, although recent work on tuned range- 
separated hybrids suggests that both quasiparticle gaps 
and outer valence spectra from a gKS approach can be 
in quantitative agreement with photoemission and GW 
calculations d ^ 

For these reasons, the BSE approach is increasingly be¬ 
ing used to predict excitation energies for molecules and 
is an alternative to TD-DET. However, there are, to date, 
no general assessments of the quality of BSE results for 
low-lying neutral excitations of isolated molecules. Al¬ 
though several benchmarks of TD-DET have been re¬ 
ported detailing its accuracy for different choices of xc 
functional relative to wavefunction-based methods,no 
such systematic effort has been undertaken for the BSE 
approach. Here, we evaluate the accuracy of BSE neu¬ 
tral excitations compared to values for 28 small organic 
molecules calculated by Thiel and coworkers with high- 
level wavefunction-based methods.This set of 28 
molecules, hereafter referred to as “Thiel’s set”, includes 
103 singlet and 63 triplet excitations, all computed with 
multiple coupled-cluster level methods.Eollowing prior 
studies with TD-DET, we benchmark to theoretical val¬ 
ues, rather than experimental data: the compared cal¬ 
culations employ the same basis set, all atomic positions 
are identical, vibrations and temperature effects are ne¬ 
glected, and there is no solvent or other environmental 
conditions. With this benchmark, we are able to provide 
a general assessment, as well as guidelines and rationale 
for the successful application of BSE to molecular sys¬ 
tems. 

The article is organized as follows: Sec. H is a general 
presentation of the BSE formalism. Sec. HI details the 
practical calculations and presents Thiel’s set. Sec. IV 
presents the BSE results for the Thiel’s set. In Sec. V 
we discuss the BSE results and analyses the causes of 
success and of failures. The study is concluded in Sec. VI. 
Hartree atomic units will be used throughout the text 
{h = e = ao = 1). 


II. BETHE-SALPETER EQUATION FOR MOLECULES 

The BSE is an equation for the two-particle Green’s 
function, more precisely for its electron-hole time¬ 
ordering. Rigorously derived from many-body perturba- 



FIG. 1. Typical BSE calculation flowchart. 


tion theory, the BSE, when using a static approximation, 
is completely analogous to TD-DET expressed through 
Casida’s equations.Here we describe the most salient 
features of a typical BSE solution in its most common 
practical implementation.^ 

Eigure 1 shows a flowchart of a BSE evaluation of neu¬ 
tral excitation energies. The calculation consists of 3 
main steps: a self-consistent mean-field DET calculation 
in the gKS scheme; a perturbative GW calculation to ob¬ 
tain the quasiparticle energies; and a final BSE solution 
to produce the excitation energies and the correspond¬ 
ing oscillator strengths /g. 

The initial DET step is denoted gKS since the starting 
mean-field might be standard Kohn-Sham with a local 
xc approximation, or it might be derived from a more 
general hybrid functional or even Hartree-Eock. This 
first step produces eigenvalues and wavefunctions that 
are used to evaluate the screened Coulomb interaction 
W and the GW self-energy. 

Details of the calculation of the GW self-energy for 
atoms and molecules can be found in Ref. 42. Let us 
simply write the self-energy E as 

S(l,2)=zG(l,2)ty(l+,2), (1) 

where the composite index 1 is short for position, time, 
and spin (ri,ti,cri). 1+ indicates the limit as time goes 
to ti from above. G is the one-particle Green’s function 
and W is the screened Coulomb interaction. The GW 
self-energy produces quasiparticle energies which are, by 
definition, the binding energies of electrons or holes in a 
system. These energies are precisely the observables mea¬ 
sured by photoemission (for occupied states) and inverse- 
photoemission (for unoccupied states). In practice, the 
GW quasiparticle energies show good agreement with 
experiment, albeit with a notable starting point depen¬ 
dence. We will discuss this further below. 

In the end, BSE is a Dyson-like equation for the so- 
called two-particle correlation function L. The full equa- 
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tion reads 

L(l,2;l',2') = Lo(l,2;l',20 

J d3d4d5d6io(l,2;4,3) 


+ 




transitions, and the two types of transitions are coupled 
through the blocks B and —B. The neglect of the cou¬ 
pling B leads to the Tamm-Dancoff approximation.^^ 
The only difference between TD-DFT and BSE lies in 
the specific expression of the matrix elements in A and 
B. In BSE, if i and j are occupied states and a and 
b are unoccupied states, these elements read, for spin- 
restricted calculations. 


where the non-interacting correlation function Lq is ex¬ 
pressed as 


Lo(1,2;1',2') = C?(1,2')G(2,1') (3) 


-a^/'^{ia\jb) + W,f{oj = 0) (9a) 

Bia = -a^^^{ia\bj) + = 0), (9b) 


and M is simply the sum of the Hartree potential and 
the self-energy 

M(3,4)=i;^(3)(5(3,4) + E(3,4). (4) 

When the indices are contracted, L and Lq yield the usual 
interacting and non-interacting polarizabilities 

X(l,2) = -zL(l,2;l+,2+) (5a) 

Xo(l,2) = -zLo(l,2;l+,2+). (5b) 

When expressed in this form, the BSE in Eq. (2) and 
the central equations of TD-DET in the linear response 
formalism 

X(l,2)=xo(l,2)+ y'd3d4xo(l,3)^^^X(4,2) (6) 

are linked in rather intuitive fashion. 

In practice, the BSE is generally solved using the 
screened Hartree-Eock approximation to E, a choice that 
can alternatively be viewed as a GW approximation to 
E in the static limit. Hence, the BSE kernel simplifies to 
the following frequency independent expression: 


+ iW(r 3 ,r 4 ,uj = 0)S(r3 - r 5 )S(r 4 - re), (7) 

where v is the bare Coulomb interaction in the previous 
equation. 

With this static assumption, the BSE can be recast 
into a matrix form in a transition space spanned by the 
orbital products (pi(r)(pj(r) where pairs of states i and j 
are either occupied/unoccupied or unoccupied/occupied. 
Thus the BSE results in an eigenvalue problem with the 
same block form as the TD-DET equations 



where (iajjb) are Coulomb integrals in Mulliken notation. 
The coefficient is set to 2 in the case of a singlet 

final state or to 0 in the case of a triplet final state. 

The eigenvalue problem posed by the BSE as shown in 
Eq. (8) is numerically cumbersome: the matrix size grows 
as the square of the number of atoms (2 times the num¬ 
ber of occupied states times the number of unoccupied 
states) and it is furthermore a non-symmetric eigenvalue 
problem. However it is well known from TD-DET that 
this problem can be reduced to a symmetric eigenvalue 
problem whose size is cut in half.^^’^^ After some algebra, 
the problem can be recast as 


cz, = ( 10 ) 


where C = {A — By/^{A + B)(A — B)^!'^ is a symmet- 
ric matrix that is half the size of the initial problem in 
Eq. (8). The above expression assumes matrix {A — B) 
to be positive definite. Erom the knowledge of an eigen¬ 
vector Zg, one can build both Xs and as 


= - 
^ 2 


{A - Bf/"^ + Q,s{A - B)-^/^ 
{A - Bf/"^ - ^s{A - B)-^/'^ 


Zs 

Zs- 


(lla) 

(llb) 


In fact, here the calculation of the square root of matrix 
{A — B) requires another diagonalization. Note that in 
TD-DET, (A—B) is a diagonal matrix and its square root 
is readily obtained. However recent work^^ has proven 
that Cholesky decompositions can be a work-around to 
avoid this second diagonalization. 

A TD-DET calculation would be essentially identical 
to this except that the eigenvalues entering in A would be 
gKS eigenvalues instead of quasiparticle energies from a 
GW approxomation, and except that the W term would 
be replaced by the xc kernel f^c (with a different index 
ordering). 

Having reviewed the BSE formalism, we will turn to 
the practical application of it to Thiel’s set.^^’^^ 


where fig are the neutral excitations and (X^, T^) are the 
eigenvectors. The complex conjugation has been dropped 
because the wavefunctions are assumed to be real-valued. 
Just as in TD-DET, the upper block A accounts for res¬ 
onant transitions from occupied to unoccupied orbitals, 
whereas the lower block — A accounts for the antiresonant 


III. THIEL’S SET: TECHNICAL ASPECTS 

In the present study we evaluate the quality of the 
BSE neutral excitation energies for Thiel’s set of 28 small 
organic molecules. 
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Unsaturated Aliphatic Hydrocarbons 



Aromatic Hydrocarbons and Heterocycles 



Aldehydes, Ketones and Amides 


Nucleobases 

FIG. 2. The 28 molecules contained in Thiel’s set. H is white, 
C is light blue, N is dark blue, and O is red. 


Our calculations are performed with the MOLGW 
code^^’^^ which is an implementation of GW and BSE 
many-body perturbation theory with Gaussian basis 
functions. MOLGW relies on an external library, LIBINT,^^ 
to evaluate the Coulomb integrals. The xc energies, 
potentials, and kernels for different starting gKS DFT 
mean-fields are obtained from the LIBXG library.The 
philosophy behind MOLGW is to prioritize accuracy and 
ease of development over computational efficiency, and 
thus MOLGW is suitable for small molecular systems. So 
MOLGW solves the random-phase approximation equa¬ 
tion [i.e. Eqs. (9a-9b) without the last term] for the 
spectral representation of VF, and thus computes the 
GW self-energy analytically. In contrast with other 
implementations,^^’"^^ we do not employ auxiliary basis 
functions to expand the 4-center Coulomb integrals, and 
so the final GW quasiparticle and BSE excitation ener¬ 
gies are exact within the selected basis set. 

Thiel’s set contains 28 organic molecules that consist 
of just four different elements (C, N, O, and H) with 
the largest molecule being naphthalene CioHg. The ge¬ 
ometries of the molecules in Thiel’s set were relaxed 
within MP2; coordinates for all molecules can be found 
in the supporting information of Ref. 38. The set com¬ 
prises tabulated reference excitation energies for 103 sin¬ 
glet and 63 triplet final states. The reference data have 
been obtained from several flavors of coupled-cluster 
theory, namely CC 2 , CCSD, and CC3,^^, and from 
complete-active-space second-order perturbation theory 
(CASPT 2 ).^^ The reference excitations of Thiel and 
coworkers are referred to as “Best Theoretical Estimates” 
(BTEs for short), and go beyond a weighted average 
of the different theoretical approaches. BTEs are the¬ 
oretical values that have been corrected with some hu¬ 
man intuition about the usual discrepancy between these 
methods and reality. Indeed, BTE values most often lie 
outside the range of the calculated values. Note that. 



FIG. 3. Basis set convergence of the first ^Biu excitation in 
ethene within GGSD (from Ref. 38), from TD-B3LYP, and 
from BSE based on B3LYP inputs. The excitation energy for 
the largest basis set (aug-cc-pVQZ) has been used as the zero 
for each theoretical approach 


consistent with the Thiel group’s subsequent TD-DET 
study,we disregarded the tabulated double excitation 
of tetrazine (C 2 N 4 H 2 ) which so far could not be captured 
by TD-DET or BSE. This explains why we refer to 103 
tabulated values instead of the 104 that appear in the 
original work. 

The original calculations performed on Thiel’s set used 
the so-called TZVP basis set of Alrichs and coworkers. 
This relatively limited basis was used so that the highly 
demanding calculations required to build the BTEs were 
feasible. Eor the sake of comparison, we employ the same 
basis set in our calculations in this work. The TZVP ba¬ 
sis contains 3 series’ of valence basis functions, but only 
one series of polarization functions {d orbitals for C, N, 
O, and p orbitals for H); it contains no diffuse functions. 
This basis set yields unconverged results as exemplified 
in Fig. 3 for the first excitation in ethene C 2 H 4 , which 
is one of the smallest molecules of the set. Although 
all methods considered here (coupled-cluster, TD-DET, 
and BSE) are clearly unconverged for the TZVP basis set 
compared to, for example, the Dunning aug-cc-VQZ ba¬ 
sis set,^^ it is demonstrated in Fig. 3 that the convergence 
rate is similar for the different approaches, justifying the 
use of a smaller basis. The deviation of the TZVP value 
from the converged value ranges from 0.35 to 0.45 eV 
across the different theoretical methods. Because the er¬ 
ror of all these methods with the TZVP basis agree within 
0.1 eV, we expect our calculations to trend across theo¬ 
retical schemes. The Thiel group has also shown that the 
conclusions drawn from the smaller TZVP basis remain 
valid with the larger aug-cc-pVTZ basis set that includes 
diffuse functions. 

With these preliminaries done, we are now ready to 
analyze the performance of BSE for the 28 selected 
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molecules of Thiel’s set. 


IV. PERFORMANCE OF THE BSE FOR THIEL’S SET 

As shown in Fig. 1, the BSE excitation energies rely 
on eigenenergies and wavefunctions from a prior self- 
consistent gKS DFT calculation. As mentioned, a strong 
dependence of the GW quasiparticle energies on the 
DFT starting point has previously been discussed in the 
literature;^^’^^’^^’^^”^^ thus it is not surprising that BSE 
excitation energies, which are built upon GW quasipar¬ 
ticle energies (as shown in Fig. 1) will also exhibit such 
a dependence. Although the influence of the DFT start¬ 
ing point was mentioned in earlier works,to date, no 
systematic quantitative study has been performed. 

Hereafter, we will assess the BSE via evaluation of their 
deviation from the reference BTEs of Thiel’s set for both 
singlet and triplet excitations. The BSE is solved us¬ 
ing GW quasiparticle energies that have been obtained 
from several xc approximations to the gKS DFT starting 
point. We have selected 4 different xc approximations 
that are reasonably representative of the popular choices 
for molecules. PBE^"^ is a pure semi-local functional with 
no exact exchange. B3LYP^^ is a hybrid functional con¬ 
taining 20 % exact exchange, whose 3 parameters have 
been adjusted to yield good thermodynamic data, and 
to this day, is one of the most widely-used functionals 
in the quantum chemistry community. BHLYP^^ is an¬ 
other hybrid functional due to Becke, but contains a sig¬ 
nificantly larger content of exact exchange, 50 %. This 
functional was identified as one of the best starting points 
for GW in our previous study.^^ Tuned CAM-B3LYP,^ 
that we label tCAM-B3LYP in the following, is a range- 
separated hybrid that has the correct full long-range ex¬ 
change {a /3 = 1). It is constructed to yield accurate 
results in TD-DFT. 

The BSE results will then be labeled BSE@PBE, 
BSE(aB3LYP, BSE^BHLYP, and BSE@tCAM-B3LYP 
respectively. We reiterate that even though it is not ex¬ 
plicitly stated in the short-hand notation, an intermedi¬ 
ate GW calculation is always performed. 


A. Singlet excitations 

Figure 4 illustrates the correlation between our com¬ 
puted BSE singlet excitation energies and the reference 
BTEs evaluated by Thiel’s group.Perfect agreement 
would be the case if all points were to he along the diag¬ 
onal line. BSE@PBE consistently yields singlet excita¬ 
tion energies that are too low: almost all data points are 
below the diagonal. BSE@B3LYP is much improved but 
still somewhat underestimates the excitation energies for 
this set. BSE^BHLYP and BSE@tCAM-BL3YP, how¬ 
ever, are in excellent agreement, with narrow scattering 
around the diagonal. The data in Fig. 4 remarkably fol¬ 
low the fit by a straight line, whose slope is very close 






FIG. 4. Correlation plots for singlet excitations between BSE 
using different starting points and BTE. A linear fit of the 
data is shown with a dashed line. 




FIG. 5. Singlet excitations mean signed error (upper panel) 
and mean absolute error (lower panel) for different schemes. 
TD-DFT based on B3LYP (from Ref. 39) is given as a com¬ 
parison. 


to unity. This means that for a given starting point the 
error is quite constant irrespective to the excitation en¬ 
ergy. 

Thus, whereas semi-local functionals like PBE are not 
suitable as a starting point for this set of small organic 
molecules, hybrid functionals do much better, and it ap¬ 
pears that a larger content of exact exchange improves 
the agreement with respect to the best theoretical esti¬ 
mates. 

In Fig. 5, we report the mean signed error (MSE) and 
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FIG. 7. Triplet excitations mean signed error (upper panel) 
and mean absolute error (lower panel) for different schemes. 
TD-DFT based on B3LYP (from Ref. 39) is given as a com¬ 
parison. 




FIG. 6. Gorrelation plots for triplet excitations between BSE 
using different starting points and BTE. A linear fit of the 
data is shown with a dashed line. 

mean absolute error (MAE) with respect to BTEs for the 
different approaches considered in this paper, and cite the 
TD-B3LYP error^^ as a reference. We select TD-B3LYP 
because it performs best for Thiel’s set among all TD- 
DFT xc functionals.^^ In fact, for the type of excitations 
considered in Thiel’s set - no charge transfer or Rydberg 
excitations - TD-B3LYP performs so admirably that we 
could not expect BSE to outperform it. As previously no¬ 
ticed, the results reported in Fig. 5 show a strong depen¬ 
dence of the BSE excitation error on the starting point. 
More precisely, BSE(RPBE underestimates all the exci¬ 
tation energies by almost 1 eV. BSE(RB3LYP also yields 
excitation energies that are too low. However, with a 
BHLYP or tCAM-BL3YP starting point, the BSE re¬ 
sults can indeed challenge the best TD-DFT excitation 
energies, yielding results with an MAE of around 0.25 eV. 

In conclusion, for singlet excitations, BSE with a prop¬ 
erly chosen starting point can be a predictive tool for sim¬ 
ple neutral excitations of small organic molecules. We 
will return to the starting point dependence further in 
Sec. V. 


analogous to TD-DFT for triplets. From the correlation 
plots shown in Fig. 6, we see that all BSE triplet exci¬ 
tations are too low, regardless of the initial gKS starting 
point. Once again, the excitation energies are well fitted 
by a straight line, however with a slope that departs from 
unity. The slope ranging from 0.88 for PBE to 0.97 for 
tCAM-B3LYP shows that the error is not perfectly con¬ 
stant across the excitation energies: the larger excitation 
energies have a greater error. As expected, BSE@PBE 
produces the poorest triplet excitation energies of all. 
Hybrid functionals with some exact exchange (BHLYP 
and tCAM-BL3YP) improve results relative to the BTE, 
but the quality of the calculated triplet excitation ener¬ 
gies is poorer than for singlets. 

The errors shown in Fig. 7 confirm that with the best 
starting point (BHLYP), our BSE calculations match 
TD-B3LYP in quality, but do not do better for Thiel’s 
set. For both TD-B3LYP and BSE@BHLYP, the error is 
systematic, with an underestimation of the triplet ener¬ 
gies by 0.4 eV. 


V. DISCUSSION 

A. A strong dependence on the starting point? 


B. Triplet excitations 

Thiel’s set also contains 63 triplet excitation energies, 
and we now briefly discuss this case. It is well doc¬ 
umented that TD-DFT can have trouble with triplet 
excitations:^^ no xc functional of TD-DFT has been able 
to predict triplet energies of the molecules in Thiel’s set 
at the level obtained for singlets. 

Unfortunately, our BSE calculations show a trend very 


As shown above, the quality of the BSE excitation 
energies is strongly affected by the gKS starting point. 
Here, we would like to discuss the sensitivity of the final 
BSE result to starting point, relative to TD-DFT. 

In Fig. 8 we represent the mean signed error for Thiel’s 
set as a function of the amount of exact exchange in the 
xc functional. The TD-DFT results are from Ref. 39, 
whereas the BSE results are those reported above. Both 
approaches show a noticeable dependence on the exact 
exchange content, and in the end, they are nearly equally 
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EXX content 



FIG. 8. Dependence of the mean signed error for the singlet 
excitation energies of the Thiel’s set with respect to the con¬ 
tent of exact exchange of the underlying xc functional. The 
xc functional is used for TD-DFT (red square symbols) or as 
the starting point of GW and BSE (blue diamond symbols). 


sensitive to the xc functional. 

Interestingly, the primary difference between the TD- 
DFT and BSE schemes lies in the amount of exact ex¬ 
change that minimizes the error. TD-DFT performs best 
with 20-25 % exact exchange as in B3LYP or PBEO.^^ 
On the contrary, for BSE, the best starting points con¬ 
tain much more exact exchange (around 50 %). These 
two very different optima can be rationalized when de¬ 
composing the origin of the errors in each of the two 
schemes as we discuss hereafter. 


B. Analysis of the origin of errors 

As summarized in Fig. 1, the BSE energies are ob¬ 
tained from a series of three calculations: a self-consistent 
gKS DFT calculation, a GW calculation of the quasipar¬ 
ticle energies, and finally an evaluation of the BSE. It 
is interesting to identify which step introduces the no¬ 
ticeable starting point dependence we highlighted above. 
Or, another way of posing the same question would be: 
which among (a) GW quasiparticle energies and (b) the 
BSE solution is the most sensitive to the gKS input? 

To address this question, we need “Best Theoretical 
Estimates” for the quasiparticle energies too. Although 
it is not possible to easily access all quasiparticle ener¬ 
gies, the HOMO and LUMO energies can be obtained via 
total energy differences with the so-called ASCF proce¬ 
dure. To obtain accurate results, we applied the ASCF 
procedure within CCSD(T), the well-known and stan¬ 
dard coupled-cluster method including single and dou¬ 
ble excitations with triples introduced perturbatively.^^ 
All of our coupled-cluster calculations are done with the 
GaussianOQ code.^^ We employ the same TZVP basis set 
used by Thiel, and by us so far in this study. Again, al¬ 
though the diffuse-orbital-less TZVP basis set is, strictly 



FIG. 9. Gorrelation plot for HOMO-LUMO gap from GW 
and from GGSD(T). 


speaking, inadequate for the LUMO wavefunctions for 
many of these molecules, we use it for consistency with 
the rest of the calculations performed in this work. The 
ASCF procedure requires three separate total-energy cal¬ 
culations for evaluation of the HOMO-LUMO gap: one 
for the neutral molecule, and additional calculations for 
the cation and anion. Note that the underlying Hartree- 
Fock self-consistent field calculations have been carefully 
checked against MOLGW, since the cation and anion cases 
can be challenging and quite often converge to local min¬ 
ima. 

The comparison between GW HOMO-LUMO gaps and 
CCSD(T) gaps is summarized in Fig. 9. The results are 
in line with a previous study of the ionization potentials 
of small molecules:^^ GW on PBE largely underestimates 
the HOMO-LUMO gap, whereas hybrid functionals with 
a large fraction of exact exchange do a much better job. 
Even GW on Hartree-Fock (GIU@HF) does a decent job 
for the gaps, as highlighted in Ref. 22. It is worth noting 
that the trends for the HOMO-LUMO gaps are the same 
as those for singlet excitations. 

Let us quantify this statement by plotting the corre¬ 
lation between the GW MSE and the BSE MSE for dif¬ 
ferent starting points in Fig. 10. Our results from BH- 
LYP and tCAM-B3LYP are the most adequate for both 
singlet and triplet quasiparticle energies: these are the 
closest to the origin of Fig. 10. Furthermore, the corre¬ 
lation between the GW and singlet BSE errors is almost 
perfect: the slope of a line fit to these data is 1.02. This 
means that the starting point dependence of BSE is in¬ 
herited entirely from the GW quasiparticle energies, and 
for Thiel’s set, the BSE singlet excitations are seen to 







FIG. 10. Correlation plot between the mean signed error in 
the GW HOMO-LUMO gap and in the BSE energies for dif¬ 
ferent gKS starting points. Both singlet (full symbols) and 
triplet (open symbols) excitations are represented. Dashed 
lines are hts of the singlet or triplet results. Horizontal and 
vertical lines mark the zero error lines. 


be constant shifts applied to the GW energies. It is un¬ 
expected that the details in BSE are so insensitive to 
the starting point. Indeed, the screened Coulomb in¬ 
teraction W used in Eqs. (9a),(9b) is obtained within 
the random-phase approximation of the underlying gKS 
DFT calculation, and therefore also varies with starting 
point. Moreover, W from PBE would be expected to 
lead to more significant screening than a W constructed 
from BHLYP, since the HOMO-LUMO gaps of these two 
approximations largely differ. But obviously, these differ¬ 
ences are too subtle at the short ranges at play in these 
small molecules to influence the BSE results. The BSE 
triplet results can be seen as slightly more correlated to 
the type of IT, but the slope of a line through these data 
(^0.9) is still very close to unity. 

Finally, from Fig. 10, we observe that the intercept 
of the line fit to the errors differs from zero for both 
singlets, 0.24 eV, and triplets, -0.36 eV. From this anal¬ 
ysis, we can quantify the magnitude of intrinsic errors 
to BSE. Indeed, even with a perfect GW approach that 
produces zero error compared to the CCSD(T) HOMO- 
LUMO gaps, the BSE singlet/triplet energies would still 
deviate from those of the best theoretical estimate by 0.2- 
0.3 eV. In addition, the singlets and triplets have opposite 
signs, indicating that the BSE singlet-triplet splitting can 
be expected to be systematically overestimated by about 
0.6 eV. Furthermore, this conclusion holds independent 
of starting point since the error lines are nearly parallel. 


It is enlightening to carry out the same analysis 
for TD-DFT. It is legitimate to do so because in a 
gKS scheme, be it exact or approximate, the non-local 
exchange-correlation operator does not have a derivative 
discontinuity.^^’^^ When combining this statement with 
the piece-wise linearity of the total energy, both the gKS 
HOMO eigenvalue should be equal to minus the IP and 
the gKS LUMO eigenvalue should be equal to minus the 
EA within the exact gKS scheme. There is no funda¬ 
mental difference between the IP and the EA. As a con¬ 
sequence, an accurate xc functional in the gKS scheme 
should yield frontier eigenvalues that compare well with 
the BTEs for the HOMO-LUMO gap. 

The TD-B3LYP results yield almost perfect singlet ex¬ 
citation energies with a MSE value of -0.08 eV. However, 
B3LYP produces HOMO and LUMO gKS eigenvalues 
that strongly deviate from the ASCF CCSD(T) refer¬ 
ence. The MSE for B3LYP HOMO-LUMO gap is as 
large as -5.1 eV for the molecules in Thiel’s set. Given 
our BSE results, this strongly suggests that its accurate 
singlet excitation energies can be ascribed to a significant 
cancellation of errors. Though B3LYP and other similar 
xc functionals provide a good estimate of singlet excita¬ 
tion energies, this agreement is not supported by a satis¬ 
factory theoretical basis. Having an xc approximation 
that yields both correct HOMO-LUMO gaps together 
with high quality neutral excitations is quite possible, as 
shown by the promising recent advances associated with 
the OT-RSH functional. 

In sum, the dependence of BSE excitation energies on 
the DFT starting point can be primarily ascribed to the 
underlying GW quasiparticle energies. The BSE step 
itself has a systematic bias towards an overestimation 
of the singlet energies and an underestimation of the 
triplet energies. Although both errors are rather small 
(0.2-0.4 eV), they lead to overestimated singlet-triplet 
splittings. However, the BSE error in the singlet-triplet 
splitting is of the same order of magnitude as the best 
TD-DFT schemes. Finally, the agreement of the TD- 
DFT approximation with the reference singlets has been 
shown to rely on a cancellation of significant errors, a can¬ 
cellation that may be less complete by excitations that 
deviate from the simple nature of those considered here, 
for example in molecules more complex than for the small 
molecules in Thiel’s set. 


VI. CONCLUSION 

In this article, we have evaluated the performance of 
BSE for singlet and triplet excitation energies of Thiel’s 
set of 28 small organic molecules. The quality of our BSE 
results is found to be sensitive to the chosen DFT gKS 
starting point. Semi-local starting points, such as PBE, 
fail significantly for the molecules in Thiel’s set. Hybrid 
functional starting points in general produce much bet¬ 
ter results. Among the hybrid functionals, those contain¬ 
ing a large contribution from exact exchange (BHLYP or 
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tCAM-BL3YP) perform best. We would advocate the 
use of such functionals (or similar ones) in future BSE 
studies. 

The dependence on the starting point can be connected 
to the different content of exact exchange in each func¬ 
tional. This sensitivity, though important, is not more 
significant than the one observed for xc functionals in 
TD-DFT. The performance of BSE for singlets is clearly 
superior to its performance for triplets. The same state¬ 
ment holds for TD-DFT. 

When analyzing the origin of the error with HOMO- 
LUMO gaps evaluated from ASCF based on CCSD(T), 
we found that the entire BSE starting point dependence 
originates with the GW quasiparticle energies. The de¬ 
tails of screened Coulomb interaction W used in the BSE 
kernel is not significant for the small molecules in Thiel’s 
set. Thus, maximizing the accuracy in quasiparticle ener¬ 
gies should minimize the error in BSE. However, residual 
errors in BSE are non-vanishing. For the best quasipar¬ 
ticle energies, we predict that singlets would be over¬ 
estimated by 0.2 eV and triplets would be underesti¬ 
mated by —0.4 eV. Future investigation of the limitations 
of the BSE for singlet-triplet splitting would be desirable. 
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